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ABSTRACT 

The capture of a compact object in a galactic nucleus by a massive black hole (MBH) , 
an extreme- mass ratio inspiral (EMRI) , is the best way to map space and time around 
it. Recent work on stellar dynamics has demonstrated that there seems to be a complot 
in phase space acting on low- eccentricity captures, since their rates decrease signifi- 
cantly by the presence of a blockade in the rate at which orbital angular momenta 
change takes place. This so-called "Schwarzschild barrier" is a result of the impact 
of relativistic precession on to the stellar potential torques, and thus it affects the 
enhancement on lower-eccentricity EMRIs that one would expect from resonant relax- 
ation. We confirm and quantify the existence of this barrier using a statistical sample 
of 2,500 direct-summation N— body simulations using both a post-Newtonian and 
also for the first time in a direct-summation code a geodesic approximation for the 
relativistic orbits. The existence of the barrier prevents low-eccentricity EMRIs from 
approaching the central MBH, but high-eccentricity EMRIs, which have been wrongly 
classified as "direct plunges" until recently, ignore the presence of the barrier, because 
they are driven by two-body relaxation. Hence, since the rates are significantly affected 
in the case of low-eccentricity EMRIs, we predict that a LISA-like observatory such 
as eLISA will predominantly detect high-eccentricity EMRIs. 

Key words: Stellar dynamics - Extreme-Mass- Ratio Inspirals - Schwarzschild barrier 
- post-Newtonian dynamics. 



1 INTRODUCTION 

Massive black holes (MBHs), with masses ranging from 
some 10 4 Mq to a few 10 9 Mq are very likely present in 
the centre of most galaxies. Measurements of the kinemat- 
ics of gas and stars in the central regions of nearby galaxies 
(see e.g. |de Zeeuw|2001| |Barth|2004| |Kormendy|2004| [rIcTT 



stone|2004 \ have provided us with compelling evidence. Our 
own Milky Way (MW) is the galaxy for which we have the 
strongest observational proof a central MBH. Data based on 
16 years of observ ations set the mass of the central SMBH 
to ~ 4 x 1O 6 M (|Eisenhauer et al.||2005| [Ghez et al.||2005[ 



2008| |Gillessen et al.|2009||Genzel et al.|2010| ) 

To interact with the central MBH, stars have to find 
themselves on "loss-cone" orbits, which are orbits elongated 



enough to have a very close-in pericentre (Frank & Rees 



1976 



|Lightman"fe Shapi rofl977| |Amaro-Seoane fc Spurzem 



2001). While main-sequence stars are tidally disrupted when 



approaching the central massive black hole (MBH) , compact 
objects (stellar black holes, neutron stars, and white dwarfs) 
slowly spiral into the MBH via the gradual loss in the form 
of gravitational radiation and are eventuall y swallowed af- 



ter some 



10 J 



orbits in the eLISA band( Amaro-Seoane 



et al.||2007| |2012a|b| |Amaro-Seoane||2012[). This is the best 
way to probe general relativity (Sopuerta 2010) and, thus, 
the factory of space and time around a massive black hole. 
We will also get additional information about the binary it- 
self: in particular the masses of the system and the spins of 
the MBH can be measured to a level of precision without 
any precedent (see Amaro-Seoane et al.|2007 and references 
therein). Besides, the detection will provide us with informa- 
tion about the distribution of dark objects in galactic nuclei 
and globular clusters. 

Producing EMRIs is more difficult than producing tidal 
disruptions of stars (e.g., Rees. 1988; Magorrian & Tremaine 

i r , r , r , | o TTi II , n r\,\ l h-, r ' oTr ■^ ll r.r.r, 1 
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1999} |Syer fe Ulmer||1999| |Wang fc Merritt||2004| ) because 
while disruptions require a single passage within a critical 
radius, an EMRI is a progressive phenomenon that is only 
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successful if the small compact object suffers a large number 
of close encounters with the central MBH. 

The requirement for an EMRI to be successful is not 
just to have a small periapsis, but that the gravitational 
radiation timescale is sufficiently shorter than 2-body relax- 
ation, which could have a significant impact on the periapsis. 
Although the basic ideas are explained in detail in the re- 



cent review of Amaro-Seoane (20121, we deem it important 



to briefly summarize in the next paragraphs the fundamen- 
tals of the ideas of two different categories of EMRIs, since 
this will be crucial to understand the main results of this 
paper. 

In a galactic nucleus without dissipation processes a star 
suffers "gravitational tugs" in the regime where the evo- 
lution is dominated by close encounters with other stars. 
These tugs, since they are driven by two-body relaxation, 
are random and originated by interacting with other stars 
that happen to have a very close position. The scattering 
rate is very similar in orbital energy E and angular mo- 
mentum L. If the star gets close to a very low L, which is 
statistically probable, then the picture changes: The rate at 
which the star changes L will be much shorter than that 
at which it changes E. If we introduce a dissipation term 
in the picture, e.g. gravitational waves, the star follows the 
same initial evolution and, at some point, our test star, the 
compact object, reaches the region in which it is on a very 
radial orbit. In this case, at every periapsis passage, the test 
star loses a significant amount of energy and, hence, the 
semi-major axis shrinks. If the process is efficient enough it 
becomes an EMRI. 

Stars in very radial orbits would then scatter faster in L 
than in E, so that they would end up as direct plunges, i.e. 
being swallowed by the MBH after an insignificant amount 
of gravitational wave bursts. Although their event rate is 
much larger, they do not allow us to get a full picture of 
the scenario as complete as compared with the slow inspi- 



ralof an EMRI (but see the recent work of Berry & Gair 
2012 on bursting sources). The threshold lies around ~ 10~ 
pc ( jHopman fc A lexander 2005J and this is what has led up 
to now to the thought that the EMRI event rate should 
be dominated by the physical phenomena happening in the 
innermost volume around the MBH, of radius ~ 10~ 2 pc. 

However, while this is strictly true for Schwarzschild 
MBHs, the situation for spinning MBHs drastically changes 



the narrative. Recently, Amaro-Seoane et al. ( 2012d I proved 



that for Kerr MBHs plunges do not plunge, but spend 
a very large number of cycles in the eLISA band. I.e. 
they are simply high-eccentricity EMRIs. The authors prove 
that the event rate of both high-eccentric and also low- 
eccentric EMRIs is enhanced by the spin as compared to 
the Schwarzschild case by an amount that depends on the 
specific eccentricity and inclination of the orbit. 

The fact that compact objects on a "plunge" orbit have 
been envisaged as uninteresting has led to an effort to un- 
derstand the phenomena that could lead to the creation of 
EMRIs in a volume of radius ~ 10~ 2 pc. In this volume, res- 
onant relaxation is very likely the most important process to 



lead compact objects to EMRI orbits (Hopman & Alexan- 



der 20061. Nonetheless, while in the absence of relativistic 



effects resonant relaxation is expected to change the angular 



momentum of the stellar BHs very efficiently, ( Merritt et al 



rect A— body experiments that introducing relativistic pre- 
cession effects are both a blessing and a curse for the inspi- 
ral event rate: RR is quenched at high eccentricities, result- 
ing in more inspirals than plunges. However, this quenching 
also means that, in total, fewer BHs will reach pericentres 
that are small enough to lead to an inspiral observable by 



gravitational wave detectors such as cLISA (Amaro-Seoane 
et al.||2012a|b[ ). In this paper we present a statistical study 
of the Schwarzschild barrier (SB) with a set of 2,500 direct- 
summation N— body simulations including relativistic cor- 
rections to study and quantify this effect. We implement 
the relativistic effects using a post-Newtonian formalism as 



Kupi et al. ( 2006 1 but also, and for the first time ever in 

c scheme. 



a direct-summation integrator, a 

This paper is organized as follows: In Section [2] we 
present the physical setup and the numerical methods used. 
In Section [3] we present general results and discuss the im- 
plications for the Schwarzschild barrier in Section [4] 



2 PHYSICAL SETTING AND NUMERICAL 
METHODS 

Recently, MAMWU estimated with a few direct-summation 
N-body simulations expanded with a statistical Monte-Carlo 
study that the traditional EMRI event rate is markedly de- 
creased by the presence of a blockade in the rate at which 
orbital angular momenta change takes place. This so-called 
Schwarzschild barrier is a result of the impact of relativistic 
precession on to the stellar potential torques. Although the 
authors find that some particles can penetrate the barrier, 
EMRIs are significantly suppressed in this scenario. 

In analogy to MAMWU, the setup we consider consists 
of a central MBH of mass M. surrounded by 50 stellar mass 
black holes (BHs) of mass m*. We fix the masses of the MBH 
and the stellar black holes at 



M. = 10 b M Q . 

m* = 50 M e . 



(1) 



Initially we distribute the stars in phase space following a 
distribution of the form N(a,e 2 ) dade 2 — Nodade 2 , with a 
the semi-major axis of the BHs and e their eccentricity. The 
semi-major axes range between 0.1 mpc < a < 10 mpc. 

This setup represents roughly a relaxed distribution of 



stellar mass objects around a MBH (Freitag et al. 20061 



2011 hereafter MAMW11) showed recently with a few di- 



We note that the event rates we obtain in this study are 
only applicable to this specific, idealized system where N = 
50 heavy, equal mass stellar BHs orbit a MBH. For more 
realistic estimates one would need to take into account other 
properties of galactic nuclei, such as the stellar background, 
a certain mass distribution and larger N. 

2.1 Timescales 

In order to have EMRI events, one needs BHs on orbits with 
pericentres of only a few gravitational radii (r 9 = GM./c 2 ). 
This requires the existence of physical mechanisms for driv- 
ing BHs from their initial orbits onto very eccentric ones. 

A purely Newtonian system has two different ways of 
exchanging angular momentum L (or eccentricity e) and or- 
bital binding energy E (or semi-major axes a). The first one 
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is by two-body scattering, or non-resonant relaxation (NR) 
(see e.g. Spitzer 1987 Binney & Tremaine 2008 Amaro- 



Seoane||2012 I. Every time two objects come close they un- 



dergo scattering, changing the momenta of the scattering 
partners. This very basic mechanism exists in all gravita- 
tionally interacting systems from compact star clusters to 
galaxies. The associated time-scale for changing the angular 
momentum L of a given particle by AL ~ L is 



t N r = 4.6Myr5 1/2 '' 



V 10 6 M e ) V50M Q / KIT, 



(2) 

where jV<i is the number of stars within a sphere of 1 mpc 
and a will denote the semi-major axis in mpc throughout the 
rest of the paper. The derivation of this equation, although 
trivial, can be found in MAMW11. 

In Newtonian systems with a central massive object 
which dominates the gravitational potential, however, this 
relaxation mechanism is usually dominated by resonant re- 
laxation (RR) ( Hopman fc Alexander|2006 1 . 

The torque of a spherical distribution of stars of mass 
M*(r < a) inside a sphere of radius r on a BH with semi- 
major axis a leads to a retrograde orbital in-plane precession 
on a time-scale 



2tt 



/M.a 3 \ 1/2 

™ ~ g(e)M*(r<aj\G~ 



(3) 



where the eccentricity dependent function g(e) is given by 



9(e) = 



1 + VT 

2yr= 



(4) 



^From Eq. Q we can derive the time-scale of changes in the 
angular momentum, which for our system is given by 



trr 



5.9 x 10 4 yr / M. \ 



©7 V50M Q 



,3/2 



(5) 



where /3 S is a factor of order unity. For the derivation of 
Eq. (5} we defer the reader to reference MAMW11. This 
mechanism for relaxation, RR, is much more efficient than 
NR at these distances because the particles interact through 
coherent torques in resonant Keplerian orbits. It therefore 
could lead in principle to an enhancement in the EMRI event 
rate. 

Relativistic effects introduce two new time-scales. The 
conservative Schwarzschild precession, appearing in partic- 
ular at the first and second post-Newtonian orders, causes 
a precession of the pericentre by an angle 



<5$ 



3nGM. 1 



(6) 



c 2 a(l — e 2 ) 

per orbit. This leads to the following associated time-scale 



27TC 2 



5$ v ' 3(GM.) 3 / 2 



5/2/1 2 

a ' (1 — e 



(7) 



where P(a) — 27r(a 3 /GM.) 1 ^ 2 is the orbital period. In a 
more convenient notation this yields 



T"SS 



(2 x 10 4 yr)a 5/2 (i- 



(8) 



M. N-3/2 

, 10 6 M e ) 

The second important time-scale is the inspiral time tgr 
via gravitational radiation only which, for high eccentricities 
(e ~ 1), is given by (Pet ers|1964 1 



^ 5c 5 a 4 _ ., 

TGW ~256G 3 m*M.(m* + M.) { C ' 



2,7/2 



(1.16 x 10 13 yr)a 4 (l 



^7/2 



/ m* \-i/ M. \ 
V50M„/ V10 6 M„/ 



10 6 Mq , 



(9) 



(10) 



for M, 3> m*. This timescale is highly sensitive to the ec- 
centricity and semi-major axis and for a typical BH in the 
system much longer than any other relevant time. However, 
for particles very close to the central MBH, gravitational 
radiation may drive them gradually into the capture radius 
leading to an "inspiral event" . 



2.2 A direct-summation code with 

post-Newtonian and geodesic corrections 

In order to integrate the initial configuration over time 
we use the publicly available planet code by Sverre 
Aarseth ( Aarseth|19 99 2003), a direct summation N— body 
integrator. We have modified this code in order to intro- 
duce relativistic corrections to the Newtonian acceleration 
I jAmaro-Seoane et al.|2012c] ) . For the studies that we present 
here we have considered the following types of dynamics: 



• purely Newtonian dynamics 

• post-Newtonian (PN) corrections 

• relativistic geodesic equations for motion of the parti- 
cles around the MBH 

In the purely Newtonian case, the integration is obvi- 
ously done without modifications to the acceleration equa- 
tions. In the second case we add the PN corrections in the 
following way: 

Newt pcriapsis shift energy loss neglected 



where the individual Fi's denote the different PN correc- 
tions to the total force on a particle, which can be found in 
Appendix [B] 

Given the high mass ratios involved in EMRIs, their 
motion around a MBH can be also approximated by solv- 
ing the geodesic equations of motion, neglecting in this way 
dissipative effects due to gravitational wave emission and 
higher-order corrections in the mass ratio. In our case, the 
geodesic equations describe the exact trajectory of a test 
mass particle around a Schwarzschild MBH. Unlike the PN 
approximation, the geodesic equations are valid even in the 
last few r g during a plunge or inspiral, however only in the 
limit m*/M. — > 0. Some orbits are expected to migrate to- 
wards plunge or inspiral orbits at pericentre distances of 
r v < 15 r g , where the errors of the PN approximation can 
already be quite significant ( Yunes & Berti 2008 1. In order to 
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test the existence of the Schwarzschild barrier at small dis- 
tances, we have implemented these corrections in the planet 
code. 

Since the geodesic equations do not contain dissipative 
terms we compare the results of using them with the conser- 
vative PN implementations, i.e. setting the dissipative cor- 



rection F 5 = in Eq. (Ill 



At any given time all the active acceleration correc- 
tions are computed only between the MBH and a stellar 
BH object. The semi-major axis and eccentricity evolution 
is tracked by monitoring the distances at periapsis r m i n and 
apoapsis r max using the standard relations 



+ ' mm 



2 

^max * min 
^max ~T~ ^*min 



(12) 
(13) 



which are valid for any acceleration correction. This calcu- 
lation does not require a PN expansion of the Keplerian ex- 
pressions for a and e and is thus consistent with the purely 
Newtonian, PN and geodesic equations of motion. 

We record a merger event whenever a particular has an 
instantaneous separation to the MBH of r < 6r g (i.e. the 
Schwarzschild last stable orbit for circular orbits). We note 
that MAMW11 use r < 8r g and thus might classify a few 
events that we call inspirals as plunges. We run Np = 500 
simulations of the models described above, with different 
random seeds for the distribution of stars, for a few Myr or 
until a merger event happens. 



3 RESULTS FROM THE SIMULATIONS 

For all the different sets P of simulations we will compute 
the average time tx,p for the occurrence of a certain event 
X in the simulation set P in the following way: 



Np 

TX,P=J2 Ti ' 
i-1 



(14) 



where the sum runs over all the total duration times T^p 
of the Np simulations and Kx,p gives the total number 
of events X in the set P. Every simulation has been given 
a burn-in time of 10 4 yr, which is of the order of trr ~ 
6xl0 4 yr (see Eq.[E) in order to discard merger events due to 
particles being created extremely close to the MBH or even 
within the capture radius by the randomization routine. 

In Eq. ( |14[ ), the possible events X can be either i for an 
inspiral event or p for a plunge event. We define the error 
for our results as the Poisson error, 



<?x,p = \J K x ,p ■ 



(15) 



The criterion for an event to be an inspiral event is taken to 
be 

(i) a cap < 1 mpc and 

(h) a cap < 1.5a(t ca p - 500 yr) . 

The second condition ensures that the semi-major axis of the 
merging body (a cap ) has shrunk significantly prior to cap- 
ture in order to dismiss plunges with low semi-major axis. 



The choice of 500 years has empirically proven to distinguish 
perfectly between plunges and inspirals. 



3.1 Newtonian Simulations (Set SI) 

In this section we present the results for our study using 
only purely Newtonian accelerations, i.e. only the term Fq in 



Eq. (111. In this case, the individual objects exchange energy 



and angular momentum efficiently via resonant relaxation 



(RR). For the plunge time, using Eq. (141, we find 



T p j = (3.7 ±0.2) x 10 4 yr. 



(16) 



This agrees very well with the RR timescale given by Eq. |5| 
for a typical particle, which confirms that this is the domi- 
nant mechanism at these radii for driving stellar objects into 
the central body in the purely Newtonian case. Of course, we 
do not identify any inspirals in the absence of gravitational 
radiation effects. 



3.2 Simulations including 2.5PN corrections (Set 
SII) 

We now add only the dissipative effects due to gravitational 
radiation emission, which appear at 2.5 PN order, to the 
acceleration equations. 

The analysis of the simulations now gives 



t p ,u = (3.8 ±0.2) x 10 4 yr, 
n,n = (2.3 ±0.4) x 10 5 yr. 



(17) 
(18) 



Now, with the inclusion of the effects of gravitational ra- 
diation, gradual inspirals into the MBH are possible. This 
converts a subset of the plunge events from the simulations 
in set SI into inspiral events. However, the inspiral time tgw 
given by Eq. ( |10[ ) only becomes smaller than the RR time 
for 



o(l — e) < 5r 9 



(19) 



which is smaller than the assumed capture radius. Thus, 
the efficient RR still drives the majority of particles into 
the capture radius before they can decouple from the stellar 
background and undergo a clean inspiral. In other words, 
the transfer of angular momentum to more eccentric orbits 
by RR is faster than the circularization by the dissipative 
2.5PN term. 



3.3 Simulations including 1PN, 2PN and 2.5PN 
corrections (Set SHI) 

In this set of simulations we include all PN terms up to 
2.5PN order in our calculations. This introduces prograde 
Schwarzschild precession in addition to the dissipation pro- 
duced by gravitational-wave emission. This effect is expected 
to increase the associated times for inspiral and plunges, 
since it eliminates efficient RR at high eccentricities. We 
find 



T Pt m = (1.3 ±0.2) x 10 6 yr, 
■n.iii = (2.0 ±0.4) x 10 6 yr- 



(20) 
(21) 
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As expected, the inspiral and plunge times are now of the or- 
der of a two-body relaxation time, Eq. Q. We also see that 
the substantial difference between plunge and inspiral times 
seen in SII vanishes, because now the relaxation time-scale 
is much higher and gravitational radiation can more easily 
decouple the star from the stellar background. Compared 
to the previous set now the 2.5PN term is able to circular- 
ize and shrink the BH orbit faster than classical relaxation 
increases eccentricity. 

3.4 Simulations including 1PN and 2PN 
corrections (Set SIV) 

In order to compare our PN results with the results using 
geodesic equations of motion, we also ran a set of simulations 
with only the conservative 1PN and 2PN terms. In this set 
we find a plunge time of 

r p ,iv = (6.6 ±0.5) x 10 5 yr. (22) 

This is consistent with the time we would obtain from SIII 
when combining plunge and inspiral events and shows again 
that the 2.5PN term does not change the important mech- 
anisms for angular momentum transfer. 

3.5 Simulations considering geodesic motion 
around the MBH (Set SV) 

In this set of simulations we investigate the system using 
Newtonian forces to describe the gravitational interactions 
between the stellar black holes and to describe the inter- 
action between the MBH and individual stellar black holes 
we use the exact solution of the motion of a test mass in 
a Schwarzschild metric (geodesic motion). This does intrin- 
sically exclude dissipative effects and therefore, the results 
of this subsection should be compared to those from the set 
SIV. In the limit of m* <C A/., the geodesic equations (in 
harmonic coordinates) expanded up to 2PN order and the 
conservative 2PN equations agree (see Appendix |A| , and 
hence they are consistent descriptions at this level of ap- 
proximation. 

In this set of simulations we obtain 

T p ,v = (7.1 ± 0.7) x 10 5 yr, (23) 

which is consistent with T Pl rv. This agreement means that 
the motion very close to the MBH is not relevant for the 
relaxation processes that drive BHs into plunge orbits. 



4 SCHWARZSCHILD BARRIER 

In order to quantify the nature of the "Schwarzschild bar- 
rier", we first plot the normalized presence density. This is 
a measure of the total time any particle spends in a certain 
logarithmic bin in (a, 1 — e), summed up over all particles and 
simulations. For comparison with the next figures, in figure 
[l]we depict the theoretical distribution that we can expect 
from a power-law of exponent 7 = 1.75. Next, we show the 
histogram for the Newtonian case, Fig. [2] (top panel) and the 
relativistic case (bottom panel). If we consider our specific 
setup, there are 3 different regions in the (a, 1 — e) plane 
where different mechanisms are efficient. In the right-most 



Cusp density for 7 = 1.75 




-3.6 



10"' 1(T 3 10 -2 

1-e 



Figure 1. Theoretical distribution of atruncated cusp of 7 = 1.75 

region, where pericentres are large, RR plays the dominant 
role. The left border of this region is roughly given by the 
appearance of the Schwarzschild precession which inhibits 
the BHs from experiencing coherent torques. Following the 
derivation in MAMW11, the time-scale for changes in angu- 
lar momentum due to an enclosed distribution of stars with 
mass m, acting as a coherent torque is 

The number of stars within a sphere delimited by the BHs 
semi-major axis, N(a) is related to the density profile p(r). 
For a general power law p(r) tx r the number of stars 
within a certain radius a becomes 

N(r <a) = iV < io 3 ~' r , (25) 

where iV<i is the number of stars within a sphere of radius 1 
mpc. The condition for the Schwarzschild barrier is that the 
relativistic precession time-scale, Eq. Q, becomes smaller 
than T co h, i.e. 

oCl-eY^^Mv^F^). (26) 

In our model with 7 = 2 and JV<i ~ 5, we obtain the relation 

N{r < a) w 5a, (27) 

and thus the barrier at 

a S B~C SB (l-e 2 )- 1/3 , (28) 

where Csb ~ 0.35 in this particular order of magnitude com- 
parison. This line is shown in blue in Fig. [2] The bottom-left 
side of it is the region where RR is inefficient. However, NR 
is still in play. The next delimiter is placed by the inspiral 
time-scale tgw. As soon as tgw < tnr, BHs decouple from 
the stellar background and inspiral gradually, driven by en- 
ergy loss through gravitational radiation. The condition for 
this, using Eq. (10} and Eq. @, yields 

1 E m-2/ M ' \ ( m * \~ 2/7 

x^^V^ 7 . (29) 
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Newtonian Dynamics 



For a = [1.5 - 2.0 mpc] integrated presence density 












no PN 

- only 2.5PN 

all PN 

onlylPN,2PN 
geodesic 



log (1-e) 



Dynamics with 1PN, 2PN and 2.5PN 




Figure 2. Integrated presence density for the Newtonian (set 
SI, top panel) and the relativistic case (set SIII, bottom panel). 
The shaded box marks the region of the slice analyzed in Fig. [3] 
The lines indicate the position of the Schwarzschild barrier with 
Csb = 0.35 (blue) and the limit for capture onto inspiral orbits 
for non-resonant relaxation (green). 



This line for gravitational capture against non-resonant re- 
laxation is shown in green in Fig. [2] The relativistic set 
(bottom panel) shows the characteristic inspiral lines, which 
cause the presence very close to the merger limit (a < 6r ff ) 
to be depleted in comparison to the Newtonian case, where 
the particles are scattered into plunge orbits onto the MBH 
instead. 

Fig. [3] shows the presence integrated over a small slice 
in semi-major axis (the shaded area in Fig. [2} for all the sim- 
ulated sets. In this plot we show the location of the barrier 
at Csb = 0.35, shown as the blue region. Since we average 
over a certain part in the semi-major axis, the lines from 
Fig. [2] now appear as regions marking the lower and upper 
boundary at the limiting semi-major axis values. 

The same decrease in presence density can be observed 
for the runs with no dissipative forces, set SIV and SV, which 
both show the same behavior towards high eccentricities. 
This confirms the validity of the PN approximation in this 



Figure 3. Integrated presence density for a semi-major axis slice 
1.5 mpc < a < 2.0 mpc. The shaded boxes represent the posi- 
tion of the Schwarzschild barrier with Csb = 0.35 (blue) and the 
threshold between the regime in which the evolution is dominated 
by dynamics and when the binary decays via gravitational loss 
(green) 



Series 


Plunge (yr x ) 


Inspiral (yr - 1 ) 


SI: Newtonian 


(2.7 ±0.2) x 10~ 5 




SII: Only 2.5 PN 


(2.6 ±0.2) x 10" 5 


(4.3 ±0.6) X 10" 6 


SIII: Full PN 


(8± 1) x 10- 7 


(5±1) x 10- 7 


SIV: 1PN, 2PN 


(1.5 ±0.1) x 10" 6 




SV: Geodesies 


(1.4 ±0.1) x 10" 6 





Table 1 . Comparison of the different event rates for the different 
scenarios studied. 



regime for our purposes. It also suggests that the system 
barely access the regime where the PN and geodesic dy- 
namics differ, i.e. where the relativistic precession predicted 
by the geodesic equations of motion is significantly different 
from the 1PN and 2PN corrections. 

The drop for set SIII as compared to IV and V is 
stronger, because of the presence of inspirals depleting the 
very high eccentricity region. 

Set SI and SII also show consistent results. This was ex- 
pected because the decoupling into inspiral orbits happens 
on length scales below the merger criterion for resonant re- 



laxation, see Eq. ( 19 1 



We summarize the event rates for our idealized system 
in Table [I] Our rates are only meant to compare different 
sets of simulations and can only be read as absolute values 
for a real galactic nucleus. 



5 DISCUSSION 



Recently, Merritt et al. ( |2011| estimated with a few direct- 
summation A^— body simulations expanded with a statisti- 
cal Monte-Carlo study that the "traditional EMRI" event is 



© 2012 RAS, MNRAS 000,[T|fTol 



The Schwarzs child barrier 7 



markedly decreased by the presence of a blockade in the rate 
at which orbital angular momenta change takes place. This 
so-called "Schwarzschild barrier" is a result of the impact 
of relativistic precession on to the stellar potential torques. 
Although the authors find that some particles can penetrate 
the barrier, EMRIs are significantly suppressed in this sce- 
nario. 

In this study we investigated the effects of relativistic 
corrections on the event rates for EMRIs compared to New- 
tonian dynamics using a PN approach as presented in the 



original work of Kupi et al. (2006) but also and for the first 



time with the implementation of geodesic equations. For this 
purpose, we ran different sets of 500 simulations each, for 
combined durations of ~ 100 Myr in order to obtain statis- 
tically solid results. We confirm the quenching of RR in the 
presence of Schwarzschild precession, i.e. the Schwarzschild 
barrier. Comparing full PN simulations (up to order 2.5) 
with the Newtonian ones, we find a ratio of the time-scales 
for the capture (combined plunge and inspiral) of 



tgr/ttm 



cwtoman 



21 ± 7 



(30) 



and for the absolute value for our setup an EMRI event rate 
of < lMyr" 1 . This value serves as an order of magnitude 
estimate for galactic centers similar to our idealized setup. 

In order to investigate the validity of the barrier at 
high eccentricities and very small pericentres, we have imple- 
mented the geodesic equations of motion around the MBH 
as an alternative to the PN corrections. We find that the re- 
sults we obtain from both methods are consistent although 
the relativistic precession they predict is significantly differ- 
ent near the last stable orbit. This suggests that the stellar 
dynamics of the systems we have studied does not access 
significantly the regime where the dynamics as described by 
the geodesic equations and the PN corrections is different. 
Therefore we deem it necessary to address the scenario to 
check whether these results are recovered when we increase 
the number of stars N. We hence plan on expanding the 
studied system for usage with TV— body codes. This will al- 
low us to investigate realistic galactic nuclei and thus make 
more precise statements about the absolute event rates to 
be expected. 

We corroborate the existence of the "Schwarzschild bar- 
rier" and its impact on the rate at which orbital angular 
momenta change takes place at these distances. It is im- 
portant to note that recently |Amaro-Seoane et al.| ( |2012d[ ) 
proved that for Kerr MBHs, direct plunges are in reality 
high-eccentricity EMRIs, and that the rate in enhanced de- 
pending on the spin and inclination. Although our simula- 
tions probe only an idealized case in which we study a pe- 
culiar, though representative distribution, the consequences 
are clear: The event rate of EMRIs for a LISA-like mis- 
sion such as eLISA ( Amaro-Seoane et al. 2012a|b I will be 
dominated by high-eccentric EMRIs, which are not "direct 
plunges" if the central MBH is spinning. On the other hand 
EMRIs produced in the region in which it was believed that 
resonant relaxation could represent a significant enhance- 
ment in the rates, will suffer a drastic cut in the rates and, 
hence, they will be only a negligible contribution. 
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APPENDIX A: GEODESIC EQUATIONS FOR A 
PARTICLE ORBITING A BLACK HOLE 

In this appendix we write down the geodesic equations of 
motion in a form that is suitable to be included in a N-body 
code that uses a Newtonian-type formulation of the equa- 
tions of motion. In order to compare results with the cases 
where PN corrections are used we write the geodesic equa- 
tions using harmonic coordinates for Schwarzschild, which 
are compatible with the harmonic gauge condition of PN 
theory. 

Since our particles represent stellar objects we need 
to consider the geodesies for massive particles (i.e. time- 
like geodesies). Given our system of spacetime coordinates 
{x»} = {t.ar*} Qi,v,>-- = 0-3; i,j,... = 1-3), a 
geodesic will be given by {x^(t)}, where r denotes the par- 
ticle's proper time. The components of the velocity vector 
are defined as 



_ rfs"(r) 
dr 

This four-velocity vector satisfies: 



-c , 



(Al) 



(A2) 



where g M „ is the Schwarzschild metric in our coordinate sys- 
tem and c denotes the speed of light. Since we are interested 
in geodesies, the velocity vector must satisfy the following 
equation of motion ( |Misner et al.|1973[ ). 

u !/ V„u M = 0, (A3) 

where V M denotes the canonical covariant derivative associ- 
ated with the spacetime metric <? M „. Expanding this equation 
we have 



o, 



(A4) 



being TZ V the Christoffel symbols associated with the space- 
time metric <? M1/ . They are given in terms of the metric by: 



1 aft — 2 ; 



dxP dx a dxP 



(A5) 



Using the splitting of time and space we can write the ve- 
locity vector as follows: 



d 



IX = M — 



dt 



(A6) 



© 2012 RAS, MNRAS 000,[T}{To] 



8 Patrick Brem, Pan Amaro-Seoane, Carlos F. Sopuerta 



where are the velocity components in the {t,x''} 

coordinate system: 



dt(r) 



gV(t) 



dr ' dr 
Therefore, on the trajectory of the particle we can write 



(A7) 



i dx l (t) dt i t 

U = r-^ -— = V U 



rv 



(A8) 



dt dr 

where v % are the spatial components of the velocity 
i dx l (t) 

and F is the general relativistic version of the special rela- 
tivistic gamma factor, which is given in terms of the com- 
ponents of the spatial velocity and the metric tensor as: 



(A9) 



9tt + 2 9u vl +9i^vi 



which, in the weak-field limit (g t 
has the usual expression: 



-c 2 , g tl 



0,9i 



v 2 > 

7? 



, 1 
(v 



(A10) 

« Stj ), 

(All) 



At this point, we can now adopt a Newtonian point of 
view by looking at the geodesic equations for the six quan- 
tities: {x z (t) , v l (t)} , that is, for the spatial coordinates and 
spatial velocity components. They can be written as: 



dx 1 
dt 

dt/ 
dt 



f" 



(A12) 
(A13) 



where the forces, f g , are actually forces per unit mass, i.e. 
accelerations, since they should not depend on the mass of 
the body (according to the equivalence principle). Moreover, 
these specific forces depend on the spacetime metric (and its 
first derivatives) and on v l . They can be written as 



+ 



jk 



j k 



(A14) 



[see expressions in 



Given initial conditions {x' ,vl} equations ( A12|A13 1 have 
a unique solution {x l (t), v l (t)} . Note that the c factor di- 
viding the forces, when going to the right-hand side of the 
equation (multiplying the Christoffel symbols) will cancel 
the c 2 factor in the denominator of r g 
Eqs. ( [A22] >-( (A27] |]. 

Since up to now the development has been quite 
general, let us now consider the case of a non-spinning 
(Schwarzschild) MBH black hole of mass M. . The metric 
components, in harmonic coordinates, can be written in the 
following form: 



9tt = 



1 + 



9u = , 



1 + — / r 

an = n ' Uj + v + 7 



(A15) 
(A16) 
(A17) 



where 



5 y x l xi 



n' = — . 
r 



GM m 



^,From here, the components of the inverse metric are: 



9 = 



1 + 



9 = 0, 



1 a - 

ir- n n + 

1 + ^ 



" 9 ( 



(A18) 

(A19) 
(A20) 
(A21) 



where x t = Sij x 3 and n i = 8 i: j n 3 . 

The important thing to determine the forces is the 
computation of the Christoffel symbols. ^From their defi- 



nition ( A5 1 we find the following result 



Tit =0, 

1 ti — 



rii 



r' — 



r 1 



1 + 



r«=o, 



r}k — 



1 - 



2n u ( 6 



1) - 



(A22) 
(A23) 

(A24) 
(A25) 

(A26) 
(A27) 



And this determines completely the geodesic equations of 



motion in Eqs. ( A12 1 and I A13 1 



Finally, we can make a post-Newtonian expansion of 
the equations of motion. That is, an expansion for r g /r <g; 
1 , and v/c <C 1 . In our case, the expression for the force 
simplifies to [see Eq. ( |A14[ ) and Eqs. ( |A22| |-( |A27] )]: 

fl = -v\ t + 2 « l rly - r> V . (A28) 

Expanding this we get: 

f i _ _ V" 2 
J a 



t f r \ 2 / r \ 3 



1-4^+9 

r \ r 



+ 2^- 



1 + 



rijir \ v 
c 



-{n (5 ]k -n 3 n k ) - 1+ 



x--, (A29) 
c c 

where the first two rows correspond to the first two terms in 



Eq. ( A28 1 . We have expanded in Taylor series the functions 
of r g /r up to order (r g /r) 4 . We can now collect the terms 



© 2012 RAS, MNRAS 000,[T]{TO 



The Schwarzs child barrier 9 



and we find the following expression, which is valid to order 
2PN [see Eq. Jm} below]: 



f 



where 



and 



GM. r 
n ■ v 



GM. 



w. r, , . , i 

^— |(-4ipn + -4 2 pn) n 



H (£ipn+#2pn)-) , 

c c J 



n ■ v x dx 1 dx 2 

c cr dt 2cr dt 

2 s ' j 

V —V ■ V = OijV V , 



1 dr r 
2cr dt c 



Alpn —4— 5- , 

r c A 



#ipn =4, 



2^. 

r 



(A30) 



(A31) 



(A32) 

(A33) 
(A34) 
(A35) 



2006), Eq. (131)] for m* / 

2 



, 3 /n-w 

-4 ipn 

2 V c 



l + 3i/)^ + (4 + 2i/)-^ 



(B2) 



■A', 



15 „ . /n • v 



+i/ (3 - 4u) 
R 



3 / n • v \ 2 v 



+ - 



25 2 \ fn-v 



87 \ i? 



.A' 



2. 5 PN-y — ^ + — ^ 



B' 1PN =4-2i/, 

S' 2 pn = - §1/ (3 + 2m) (^) 2 + v + 2, 



2+ — +4, 



2 \ -^g 



B' 2 



24 ( R 



5 \ r 



8 -RgV^ 

5 r c 2 



(B3) 

(B4) 
(B5) 

(B6) 
(B7) 



where v is the symmetric mass ratio, v = m*M./M 2 , 
and R g = GM/c 2 . One can verify that the coefficients in 



Eq. (|A32J to Eq. (|A35j) agree with Eq. (|B2] to (|B7j for 

K = 0. 
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